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Abstract 

This paper proposes the Ricci-flow equation as a general geometric 
model for various reaction-diffusion systems (and related dissipative soli- 
tons) in mathematical biology. More precisely, we propose a hypothesis 
that any kind of reaction-diffusion processes in biology, chemistry and 
physics can be modelled by the combined geometric-diffusion system. In 
order to demonstrate the validity of this hypothesis, we review a number 
of popular reaction-diffusion systems and try to show that they can all 
be subsumed by the presented geometric framework of the Ricci flow. 

Keywords: Ricci flow, bio-reaction-diffusion, dissipative solitons and 
breathers 



1 Introduction 

Parabolic reaction-diffusion systems are abundant in mathematical biology. 
They are mathematical models that describe how the concentration of one or 
more substances distributed in space changes under the influence of two pro- 
cesses: local chemical reactions in which the substances are converted into each 
other, and diffusion which causes the substances to spread out in space. More 
formally, they are expressed as semi-linear parabolic partial differential equa- 
tions (PDEs, see e.g., [SS]). The evolution of the state vector u(x, t) describing 
the concentration of the different reagents is determined by anisotropic diffusion 
as well as local reactions: 

a t u = DAu + R(u), (d t =d/dt), (I) 

where each component of the state vector u(x, t) represents the concentration of 
one substance, A is the standard Laplacian operator, D is a symmetric positive- 
definite matrix of diffusion coefficients (which are proportional to the velocity of 
the diffusing particles) and R(u) accounts for all local reactions. The solutions 
of reaction-diffusion equations display a wide range of behaviors, including the 
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formation of travelling waves and other self-organized patterns like dissipative 
solitons (DSs). 

On the other hand, the Ricci flow equation (or, the parabolic Einstein equa- 
tion), introduced by R. Hamilton in 1982 |16| . is the nonlinear heat-like evolu- 
tion equation^ 

d t gij = -2Ri h (2) 

for a time-dependent Riemannian metric g — gij (t) on a smooth rea0 n— manifold 
M with the Ricci curvature tensor Rij This equation roughly says that we can 
deform any metric on a 2-surface or n— manifold by the negative of its curva- 
ture; after normalization (see Figured]), the final state of such deformation will 
be a metric with constant curvature. However, this is not true in general since, 
in addition to the presence of singularities, the limits could be Ricci solitons 
(see below). The factor of 2 in (|2|) is more or less arbitrary, but the negative 
sign is essential to insure a kind of global volume exponential decay^ since the 
Ricci flow equation |2]) is a kind of nonlinear geometric generalization of the 

1 A current hot topic in geometric topology is the Ricci flow, a Riemannian evolution 
machinery that recently allowed G. Perelman to prove the celebrated Poincare Conjecture, a 
century— old mathematics problem (and one of the seven Millennium Prize Problems of the 
Clay Mathematics Institute) - and won him the 2006 Fields Medal (which he declined in a 
public controversy) |41| . The Poincare Conjecture can roughly be put as a question: Is a closed 
3-manifold M topologically a sphere if every closed curve in M can be shrunk continuously 
to a point? In other words, Poincare conjectured: A simply-connected compact 3-manifold is 
diffeomorphic to the 3— sphere S 3 (see e.g., |63p. 

2 For the related Kahler— Ricci flow on complex manifolds, see e.g., |27l 128] . 

3 This particular PDE J2j was chosen by Hamilton for much the same reason that A. 
Einstein introduced the Ricci tensor into his gravitation field equation, 

1 

where Tjj is the energy— momentum tensor. Einstein needed a symmetric 2-index tensor which 
arises naturally from the metric tensor and its first and second partial derivatives. The 
Ricci tensor Rij is essentially the only possibility. In gravitation theory and cosmology, the 
Ricci tensor has the volume-decreasing effect (i.e., convergence of neighboring geodesies, see 

ED)- 

4 This complex geometric process is globally similar to a generic exponential decay ODE: 

x = -Xf(x), 

for a positive function f(x). We can get some insight into its solution from the simple expo- 
nential decay ODE, 

x = ~Xx with the solution x(t) = x$e~ At , 

(where x = x(t) is the observed quantity with its initial value xq and A is a positive decay 
constant), as well as the corresponding nth order rate equation (where n > 1 is an integer), 

x = —Xx n with the solution = + (n — 1) Xt. 

7,71— 1 ~>~n — l 
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standard linear heat equatior^ 

d t u = Au. (3) 

Like the heat equation the Ricci flow equation @ is well behaved in forward 
time and acts as a kind of smoothing operator (but is usually impossible to solve 
in backward time). If some parts of a solid object are hot and others are cold, 
then, under the heat equation, heat will flow from hot to cold, so that the object 
gradually attains a uniform temperature. To some extent the Ricci flow behaves 
similarly, so that the Ricci curvature 'tries' to become more uniform [33], thus 
resembling a monotonic entropy growth^ d t S > 0, which is due to the positive 
definitencss of the metric g^ > 0, and naturally implying the arrow of time 

|EII|3S1|2Z|. 

In a suitable local coordinate system, the Ricci flow equation ([2]) has a non- 
linear heat-type form, as follows. At any time t, we can choose local harmonic 
coordinates so that the coordinate functions are locally defined harmonic func- 
tions in the metric g(t). Then the Ricci flow takes the general form (see e.g., 

H) 

d t 9ij = Am5« + Qij(9, 9g), (4) 

where Am is the Laplace-Beltrami operator ([5]) and Q = Qij(g, dg) is a lower- 
order term quadratic in g and its first order partial derivatives dg. From 
the analysis of nonlinear heat PDEs, one obtains existence and uniqueness of 
forward-time solutions to the Ricci flow on some time interval, starting at any 
smooth initial metric go . 

The quadratic Ricci flow equation ([4]) is our geometric framework for gen- 
eral bio-reaction-diffusion systems, so that the spatio-temporal PDE (fT]) cor- 
responds to the quadratic Ricci flow PDE 

d t u = DAu + R(u) 

I I I 

dtdij = Ajuftj + Qij(g,dg) 

with: 

• the metric g = g^ on an n— manifold M corresponding to the n— dimensional 
(or n— component, or ?i— phase) concentration u(x, t); 

5 More precisely, the negative sign is to make the equation parabolic so that there is a 
theory of existence and uniqueness. Otherwise the equation would be backwards parabolic 
and not have any theory of existence, uniqueness, etc. 

6 Note that two different kinds of entropy functional have been introduced into the theory of 
the Ricci flow, both motivated by concepts of entropy in thermodynamics, statistical mechanics 
and information theory. One is Hamilton's entropy, the other is Perelman's entropy. While in 
Hamilton's entropy, the scalar curvature R of the metric gij is viewed as the leading quantity 
of the system and plays the role of a probability density, in Perelman's entropy the leading 
quantity describing the system is the metric gij itself. Hamilton established the monotonicity 
of his entropy along the volume- normalized Ricci flow on the 2-sphere S 2 |18l . Perelman 
established the monotonicity of his entropy along the Ricci flow in all dimensions |51| . 
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• the Laplace-Beltrami differential operator Am, as denned on C 2 — functions 
on an n— manifold M, with respect to the Riemannian metric gij, by 



v /det(.g) dx 



1 d 



i(vS5)p«^j) 



(5) 



- corresponding to the n— dimensional bio-diffusion term DAu; and 

• the quadratic n— dimensional Ricci-term, Q = Qij(g, dg), corresponding 
to the n— dimensional bio-reaction term, R(u). 

As a simple example of the Ricci flow equations ([2|)-(j4]), consider a round 
spherical boundary S 2 of the 3-ball radius r. The metric tensor on S 2 takes the 
form 



is independent of r. The Ricci flow equation on S 2 reduces to 

r 2 = -2(n-l), with the solution r 2 (i) = r 2 (0) - 2(n - l)t. 

Thus the boundary sphere S 2 collapses to a point in finite time (see [4"5]). 

More generally, the geometrization conjecture [60] holds for any 3-manifold 
M (see below). Suppose that we start with a compact initial 3-manifold Mo 
whose Ricci tensor i£y is everywhere positive definite. Then, as Mq shrinks 
to a point under the Ricci flow ([2]), it becomes rounder and rounder. If we 
rescale the metric g^ on Mo so that the volume of Mo remains constant, then 
Mo converges towards another compact 3-manifold M\ of constant positive 
curvature (see [T6]1. 

In case of even more general 3— manifolds (outside the class of positive Ricci 
curvature metrics), the situation is much more complicated, as various singu- 
larities may arise. One way in which singularities may arise during the Ricci 
flow is that a spherical boundary S 2 — dM of an 3— manifold M may collapse 
to a point in finite time. Such collapses can be eliminated by performing a kind 
of 'geometric surgery' on the 3 manifold M , that is a sophisticated sequence 
of cutting and pasting without accumulation of time errorsj (see [52]). After a 
finite number of such surgeries, each component either: (i) converges towards a 
3-manifold of constant positive Ricci curvature which shrinks to a point in finite 
time, or possibly (ii) converges towards an S 2 x S 1 which shrinks to a circle S 1 

7 Hamilton's idea was to perform surgery to cut off the singularities and continue his flow 
after the surgery. If the flow develops singularities again, one repeats the process of performing 
surgery and continuing the flow. If one can prove there are only a finite number of surgeries 
in any finite time interval, and if the long-time behavior of solutions of the Ricci flow J2} with 
surgery is well understood, then one would be able to recognize the topological structure of 
the initial manifold. Thus Hamilton's program, when carried out successfully, would lead to 
a proof of the Poincare Conjecture and Thurston's Geometrization Conjecture 63 . 



9ij — r 9iji 

where g^ is the metric for a unit sphere, while the Ricci tensor 



Rij ={n- l)g. 
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in finite time, or (iii) admits a 'thin-thick' decomposition of [3U]. Therefore, 
one can choose the surgery parameters so that there is a well defined Ricci flow 
with surgery, that exists for all time [52] . 

In this paper we use the evolving n— dimensional geometric machinery of the 
volume-decaying and entropy-growing Ricci flow g(t), given by equations (J5J)- 
P]). for modelling various biological reaction-diffusion systems and dissipative 
solitons, defined by special cases of the general spatio-temporal model fl|. 

2 Bio— reaction— diffusion systems 

In case of ideal mixtures, the driving force for the general diffusion DAu fl} 
is the concentration gradient — Vu, or the gradient of the chemical potential 
— Vzii of each species u.j, (i — 1, ...,n), giving the diffusion flux by the First 
Fick's law, 

J = -DVu. (6) 

Assuming the diffusion coefficients D to be a constant, the Second Fick's law 
gives the linear parabolic heat equation, 

d t u = DAu, (7) 

while, in case of variable diffusion coefficients D, we get (slightly) more general 
parabolic diffusion equation, 

8 t u = V- (DVu) , (8) 
which is still analogous to the 'linear' part of the quadratic Ricci flow equation 

dtgij = AMgij, 

due to general 'diffusion properties' of the Laplace-Beltrami operator Am- 

The n— dimensional diffusion coefficient D = D(T) at different temperatures 
T can be approximated by the Arrhcnius exponential-decay relation, 

D(f)=D„e-^, 

where Do is the maximum possible diffusion coefficient (at infinite temperature 
T), Ea is the activation energy for diffusion (i.e., the energy that must be 
overcome in order for a chemical reaction to occur) and r is the gas constant. 

Using the First Fick's first law © , the diffusion equation ([5]) can be derived 
in a straightforward way from the continuity equation, which states that a change 
in density in any part of the system is due to inflow and outflow of material 
into and out of that part of the system (effectively, no material is created or 
destroyed), 

9 t u + V-j = 0, 
where j is the flux of the diffusing material. 
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The most important special case of ([7]) is at a steady state, when the con- 
centrations u do not change in time, giving the Laplace's equation, 

Au = 0, or Aui = 0, (9) 

for harmonic functions u = {ui}. 

The stochastic version of the deterministic heat equation (J7J), connected with 
the study of Brownian motion[f| is the Fokker-Planck equation (see e.g., |31)V 

8J = -d x , [D\{x l )f] + x , x3 [Dfjix*)/] , (10) 

= cS - ' ^'i 3 = dx* dx )' Wjli ere where Dj is the drift vector and Dfj the 
diffusion tensor (which results from the presence of the stochastic force). The 
Fokker-Planck equation (TTUJ) is used for computing the probability densities of 
stochastic differential equations^] 

Also, notice that the real-valued heat equation ([7]) is formally similar to the 
complex- valued Schrddinger equation (see e.g., |32j). 

ft* = (ID 

where ip — V'C 35 -'^) is the wave-function of the particle, i = and h is 

Planck's constant divided by 2ir. 

In the remainder of this section, we will review a number of particular bio- 
reaction-diffusion systems, which are likely to be subsumed by the quadratic 
Ricci flow model |@J. 

2.1 1— component systems 

2.1.1 Kolmogorov Petrovsky Piscounov equation 

The simplest bio-re-action-diffusion PDE concerning the concentration u = 
u(x, t) of a single substance in one spatial dimension, 

d t u = Ddlu + R[u), (12) 

8 Brownian motion is the random movement of particles suspended in a liquid or gas or the 
mathematical model used to describe such random movements, often called a particle theory. 
The infinitesimal generator (and hence characteristic operator) of a Brownian motion on R" 
is ^ A, where A is the Laplacian on R n . More generally, a Brownian motion on an n— manifold 
M is given by one-half of the Laplace-Beltrami operator Am ©. 

"Consider the Ito stochastic differential equation, 

dXt = (J,(Xt,t) dt + <r(Xi, t) dWt, 

where X t £ R n is the state of an n— dimensional stochastic system at time t and Wt S K m 
is the standard mD Wiener process. If the initial distribution is Xo ~ /(x, 0), then the 
probability density of the state is given by the Fokker— Planck equation dlOI I with the drift and 
diffusion terms, 

£>i(x,t) =/ij(x,t) and (x, t) = - ^ a ik (x, t)o-^(x, t). 
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is also referred to as the Kolmogorov-Petrovsky-Piscounov (KPP) equation. If 
the reaction term vanishes, then the equation represents a pure diffusion process 
described by the heat equation. In particular, the choice 

R(u) =u(l-u) 

yields Fisher's equation that was originally used to describe the spreading of 
biological populations F°l 

The one-component KPP equation (fl2|) can also be written in the variational 
(gradient) form 

= (13) 

and therefore describes a permanent decrease (a kind of exponential decay) of 
the system's free energy functional 



F 



dx. 



where V(u) is the potential such that 

(14) 

2.1.2 Swift Hohenberg equation 

The Swift-Hohenberg (SH) equation, no-ted for its pattern-forming behavior, 
is the decaying reaction-diffusion PDE, 

<9 t u = -(l + A) 2 w + i?(u), (15) 

given by the variational (gradient) equation (|13[) with the free energy functional 



F = 



V(u) + l((l + A)uf 



dxdy, 



where R(u) is given by (|14[) . while 51 is a 2 -dimensional region in which (bio)chemical 
pattern formation occurs. 



10 In addition, the effects of convection and quenched spatial disorder on the evolution of 
a population density are described by a generalization of the Fisher/KPP equation given by 



m 



dtu = D\7 2 u + Uu — qu 



2 



where u = (x, i) represents the population density, D is a spatially homogenous diffusion 
constant, U = £/(x) is a spatially inhomogeneous growth term, and q = bio d is a competition 
term (b is a competition rate and £q is the microscopic length scale at which two particles will 
compete with one another) . One simple form of inhomogeneity considered in these works is 
a 'square well' potential J7(x) which consists of a uniform space with negative growth rate 
(termed the 'desert'), in which a single region of positive growth rate (an 'oasis') is placed. 
This model has proven to be applicable to experiments with bacteria populations in adverse 
environments 45 . 
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The time derivative of the free energy F is given by 

~dV(u 



dtF = 



du 



(1 + A)u 



d t u dxdy, 



and, since the expression in square brackets is equal to the negative right-hand 
side of (fT5l) . we have 

(dtu) 2 dxdy < 0. 

Therefore, the free energy F is the Lyapunov functional that may only decrease 
as it evolves along its trajectory in some phase space. If F has no minima, 
then when the horizontal scale of the liquid container is large compared to the 
instability wavelength, a propagating front will be observed (e.g., in chemically 
reacting flames). In this case, F will decrease continuously until the front ap- 
proaches the boundary of the medium. An alternative possibility is realized 
when F has one or several minima, each corresponding to a local equilibrium 
state in time. In this case the so-called multi-stability is possible. Therefore, 
the limit behavior of gradient systems of the form of (fTi)) is characterized by 
either a steady attractor or propagating fronts [55] , 



2.1.3 Ginzburg Landau equation 

One of the most popular models in the pattern-formation theory is the complex 
Ginzburg-Landau equation (see e.g., |56j). 

d t A = eA + (1 + ia)AA - (1 + i(3)\A\ 2 A, (16) 

where A is the complex wave amplitude, i = s is the super-criticality 

parameter, while a and /3 measure linear and nonlinear dispersion (the depen- 
dence of the frequency of the waves on the wave- number), respectively. The 
equation (|16p describes a vast array of phenomena including nonlinear waves, 
second-order phase transitions, Rayleigh-Benard convection and superconduc- 
tivity. The equation describes the evolution of amplitudes of unstable modes for 
any process exhibiting a Hopf bifurcation, for which a continuous spectrum of 
unstable wave-numbers is taken into account. It can be viewed as a highly gen- 
eral normal form for a large class of bifurcations and nonlinear wave phenomena 
in spatially extended systems 

1 The extension of the complex Ginzburg— Landau equation JTfJ, which describes strongly 
resonant multi-frequency forcing of the form 

F = fie*"* + f 2 e 2i " i + f 3 e 3iw * + c.c. 

was recently proposed in j9] by considering the analogous center— manifold reduction of the 
extended dynamical system in which the forcing amplitudes fi, f 2 , and fg are considered as 
dynamical variables that vary on the slow time scale t. Under time translations T T : A — > 
Ae*" T , they transform as /i — > /ie l " T , f 2 — > /2e 2lt ^ T , fa — * /3e 3 "* ;T . To cubic order in A the 
most general equation that is equivariant under T T is then given by 

8 t A = ai + a 2 A + a 3 AA + a 4 A\A\ 2 + a 5 A + a 6 A 2 , (17) 
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In particular, if we put a = f3 = in (HHJ), we get the real, or dissipative, 
Ginzburg-Landau equation, 

d t A = eA + AA - \A\ 2 A, (18) 

which is a gradient equation: dtA = —5F/5A, with the free energy functional 



F = - 



Using the fact that 



£ \ A \2_ l|^|4 + (Vj4) 2 



F = — \d t A\ 2 dxdy < 0, 



dxdy. 



solutions of (fT5|) at t — ► oo are either stationary field-distributions satisfying, 
for e = 1, the equation 

AA + A- \A\ 2 A = 0, (19) 

of fronts whose propagation is accompanied by a decrease of the functional F. 
The functional must reach its minimum at stable stationary solutions of (Q2 



2.1.4 Neural field theory 

The dynamical system from which the temporal evolution of neural activation 
fields is generated is constrained by the postulate that localized peaks of activa- 
tion are stable objects, or formally, fixed-point attractors. Such a field dynamics 
has the generic form [58] 

rdtu = —it + resting level + input + interaction, (20) 

where u = u(x, t) is the activation field defined over the metric dimension x and 
time t. The first three terms define an input driven regime, in which attractor 
solutions have the form 

u(x,t) = resting level + input. 

The rate of relaxation is determined by the time scale parameter r. The inter- 
action stabilizes localized peaks of activation against decay by local excitatory 
interaction and against diffusion by global inhibitory interaction. In Amari's 
formulation [3] the conceptual model (|20[) is specified as a continuous model for 
neural activity in cortical structures, 

Td t u(x,t) = -u(x,t) + h + S(x,t)+ dx'w(x-x')a(u(x',t)), (21) 



where ai = bn/i + 612/2/3, a-2 = 621 + &22 1 ^3 1 2 , as = 651/2, «6 = &6l/3. The forcing terms 
fi satisfy decoupled evolution equations on their own. In the simplest case this evolution 
expresses a de-tuning uj of the forcing fj from the respective resonance and the fj satisfy 

/; ''',/;• (7 = 1... 3). 

In general, the de-tuning introduces time dependence into i!7l . 
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where h < is a constant resting level, S(x, t) is spatially and temporally 
variable input function, w(x) is an interaction kernel and cr(u) is a sigmoidal 
nonlinear threshold function. The interaction term collects input from all those 
field sites x' at which activation is sufficiently large. The interaction kernel 
determines if inputs from those sites are positive, driving up activation (excita- 
tory), or negative, driving down activation (inhibitory). Excitatory input from 
nearby location and inhibitory input from all field locations generically stabi- 
lizes localized peaks of activation. For this class of dynamics, detailed analytical 
results provide a framework for the inverse dynamics task facing the modeler, 
determining a dynamical system that has the appropriate attractor solutions 




2.2 2— component systems 

Two-component systems allow for a much larger range of possible phenomena 
than their one-component counterparts. An important idea that was first pro- 

12 Recently, a neural attractor dynamics (NAD) was designed (see 58 ) based on a discretiza- 
tion for single neurons of Amari's neural field equation (12 1 [ i. The so— called discrete Amari 
equation describes the temporal evolution of the activity of all single neurons considering pos- 
itive and negative contributions from external input and internal neural interactions. Since 
only activated neurons can have an impact on other neurons, the neural attractor dynamics 
is nonlinear, and effects of bi-stability and hysteresis can be used for low— level memory and 
neural competition. The NAD describes the temporal rate of change of the dynamical variable 
Ui of neural activity for all behavioral neurons i. It is formulated as the following differential 
equation: 

TUl = - Ut + h + s^ h + c mot ■ o-irrn) + a^ exCjj + og^ - o&k, (22) 

where the system parameters have the following meaning: 

t, the constant relaxation rate, i.e., the time scale on which the dynamics reacts to changes; 
h, the constant negative resting level of neural activation; 

<t(.), a sigmoidal function, which maps the value of neural activity onto [0, 1], given by 
cr(u) = 1+ ^-/j u i where /3 (=100) parameterizes the slope of the resulting function; 

s^? eh , the adequate stimulus provided by sensory input of a certain duration; 
Ui, activity of behavioral neuron i, i.e., activity of behavior i; 
Cmot, a constant for weighting the motivational contribution, c mo t < \h\; 
a sclfcxc i excitatory contribution of neuron i's own activity Ui] 
a exci' au excitatory contribution of active neurons connected to neuron i; 
a inh ' 1 a " inhibitory contribution of active neurons connected to neuron i 
mi, activity of motivational neuron i, i.e., motivation of behavior i is in 1581 defined by the 
following NAD— equation, similar to (122 D : 

• . i L i „mot I „ mot I „ mot „ mot 

rim = -mi + h + Si + a sclfcxc i + a exc i - a inM , 

where 

a se°fcxc '' exc it a t° r y contribution of neuron i's own motivation m;; 

a cx° t i' au excitatory contribution of motivation neurons connected to neuron i; 

a hSii > a " inhibitory contribution of motivation neurons connected to neuron i. 

In this framework, a nonlinear neural dynamical and control system generates the tem- 
poral evolution of behavioral variables, such that desired behaviors are fixed-point attractor 
solutions while un-desired behaviors are repellers. 

This kind of attractor & repeller dynamics | 26 | provides the basis for understanding cogni- 
tion, both natural and artificial 33 29l 130] . 
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posed by A. Turing is that a state that is stable in the local system should 
become unstable in the presence of diffusion [61]. This idea seems counter- 
intuitive at first glance as diffusion is commonly associated with a stabilizing 
effect. However, the linear stability analysis shows that when linearizing the 
general two-component system 

d t u\_(D u \ / d xx u \ / F(u, v) 
d t v J \ D v J \ d xx v J + \ G{u, v) 

and perturbing the system against plane waves 

- f +\ ( "(*) 
Uk(x,*)=( - t 



close to a stationary homogeneous solution one finds [30j 

d t u u {t) \ _ ,2 / D u u k (t) \ , ( fi k (t) 
$«k(t) )~ V A,«k(t) J \ «k(*) 

Turing's idea can only be realized in four equivalence classes of systems charac- 
terized by the signs of the Jacobian R' of the reaction function. In particular, 
if a finite wave vector k is supposed to be the most unstable one, the Jacobian 
must have the signs 



This class of systems is named activator— inhibitor system after its first repre- 
sentative: close to the ground state, one component stimulates the production of 
both components while the other one inhibits their growth. Its most prominent 
representative is the FitzHugh-Nagumo equation (flZij]) . 

2.2.1 Brusselator 

Classical model of an autocatalytic chemical reaction is Prigogine's Brusselator 
(see e.g., [54]) 

d t u = DlAu + a + u 2 v - (1 + f3)u, d t v = D 2 v Av - u 2 v + 0u, (23) 

which describe the spatio-temporal evolution of the intermediate components 
u and v, with diffusion coefficients D u and D V1 while reactions 

a — ► u, Zu + v — > 6u, p + u — > v + a, u — ► c 

describe the concentration of the original substances a and /?, for which the final 
products c and d are constant when all four reaction rates n equal unity. 
A discretized (temporal only) version of the Brusselator PDE (|2"3"|) reads 

u = a + u 2 v — (1 + 0)u, v = j3u — u 2 v. 

The Brusselator displays oscillatory behavior in the species u and v when 
reverse reactions are neglected and the concentrations of a and (3 are kept 
constant. 
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2.2.2 2 component model of excitable media 

Turbulence of scroll waves is a kind of spatio-temporal chaos that exists in 
3-dimensional excitable media. Cardiac tissue and the Belousov-Zhabotinsky 
reaction are examples of such media. In cardiac tissue, chaotic behavior is 
believed to underlie fibrillation which, without intervention, precedes cardiac 
death. Fast computer-simulation of waves in excitable media have been often 
performed using the 2-component Barkley model of excitable media [5] , 



where e is a small parameter e <C 1 characterising mutual time scales of the 
fast u and slow v variables, and a and b specify the kinetic properties of the 
system. Parameter b determines the excitation threshold and thus controls the 
excitability of the medium. The term h(t) represents an 'extra transmembrane 
current'. 

Suppression of the turbulence using stimulation of two different types, 'mod- 
ulation of excitability' and 'extra transmembrane current' was performed in [46] . 
using the Barkley model (j2~4")l . With cardiac defibrillation in mind, the authors 
used a single pulse as well as repetitive extra current with both constant and 
feedback controlled frequency They show that turbulence can be terminated 
using either a resonant modulation of excitability or a resonant extra current. 
The turbulence is terminated with much higher probability using a resonant 
frequency perturbation than a non-resonant one. Suppression of the turbulence 
using a resonant frequency is up to fifty times faster than using a non-resonant 
frequency, in both the modulation of excitability and the extra current modes. 
They also demonstrate that resonant perturbation requires strength one order 
of magnitude lower than that of a single pulse, which is currently used in clinical 
practice to terminate cardiac fibrillation. 

2.2.3 Gierer Meinhardt activator inhibitor system 

Spontaneous pattern formation in initially almost homogeneous systems is com- 
mon in both organic and inorganic systems. The Gierer-Meinhardt model [14) 
is a reaction-diffusion system of the activator-inhibitor type that appears to 
account for many important types of pattern formation and morphogenesis ob- 
served in biology, chemistry and physics. The model describes the concentration 
of a short-range autocatalytic substance, the activator, that regulates the pro- 
duction of its long-range antagonist, the inhibitor. It is given as a 2-component 
nonlinear PDE system, 

dta = —p a a + pa 2 /h + D a d x ia + p a , dth = —p,hh + pa 2 + D^d^h + ph, 

where a is a short-range autocatalytic substance, i.e., activator, and h is its 
long-range antagonist, i.e., inhibitor, dta and dth describe respectively the 
changes of activator and inhibitor concentrations per second, p a and ph are 





(24) 
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the corresponding decay rates, while D a and are the corresponding diffu- 
sion coefficients, p is a positive constant. p a is a small activator-independent 
production rate of the activator and is required to initiate the activator auto- 
catalysis at very low activator concentration, e.g., in the case of regeneration. 
A low baseline production of the inhibitor, /?/,, leads to a stable non-patterned 
steady state; the system can be asleep until an external trigger occurs by an 
elevation of the activator concentration above a threshold 42J. 

2.2.4 Fitzhugh Nagumo activator inhibitor system 

An important example of bio-reaction-diffusion systems, frequently used in neu- 
rodynamics, is the 2-component Fitzhugh-Nagumo activator-inhibitor system 
[HUT] (see also 

r u d t u = D 2 u Au + f{u) - av, T v d t v = D 2 v Av + u-v, (25) 

with f(u) = Xu — u 3 ~k, which describes how an action potential travels through 
a nerve, D u and D v are diffusion coefficients, r u and t v are time characteristics, 
while k, a and A are positive constants. In matrix form, system (|25p reads 

/ T u d t u \_(Dl \ f Au \ ( \u - u 3 - k - av \ 
\ Tv d t v ) \ D 2 V )\ Av ) + \ u-v J' 

When an activator-inhibitor system undergoes a change of parameters, one 
may pass from conditions under which a homogeneous ground state is stable to 
conditions under which it is linearly unstable. The corresponding bifurcation 
may be either a Hopf bifurcation to a globally oscillating homogeneous state with 
a dominant wave number k — or a Turing bifurcation to a globally patterned 
state with a dominant finite wave number. The latter in two spatial dimensions 
typically leads to stripe or hexagonal patterns. 

In particular, for the Fitzhugh-Nagumo system (|25p . the neutral stability 
curves marking the boundary of the linearly stable region for the Turing and 
Hopf bifurcation are given by 

q*(k): l + (di + ±d 2 v )k 2 =f'(u h ), 
1n(k)-- TT%W+d 2 u k 2 =f'{u h ). 

If the bifurcation is subcritical, often localized structures (i.e., dissipative soli- 
tons) can be observed in the hysteretic region where the pattern coexists with 
the ground state. Other frequently encountered structures comprise pulse trains, 
spiral waves and target patterns. 

The reduced (temporal) non-dimensional Fitzhugh-Nagumo equations read: 

v = v(a — v)(v — 1) — w + I ai (26) 
w = bv — 7W, (27) 

where < a < 1 is essentially the threshold value, b and 7 are positive constants 
and I a is the applied current. The drift field for this model is given by 

U\{v, w) = v{a — v)(v — 1) — w, U2(v, w) = bv — jw. 
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As can be seen from (|2~T1) the null cline of the deterministic dynamics of this 
equations is the line v = jw. By substitution on the r.h.s of equation (|26[) we 
find the following equation for steady states: v(a — v)(v — 1) — — 0. 

When this system is in a noisy environment, in the limit of weak noise, 
we can approximate the dynamics of the fluctuations by the Langevin equation 

ii = v(a — v)(v — 1) v + 

7 

that is, the fluctuations run along the line v = jw. 

In particular, parameters in the FitzHugh-Nagumo neuron model [331 130) 

v = a + bv + cv 2 + dv 3 — u, u — e(ev — u), 

can be tuned so that the model describes spiking dynamics of many resonator 
neurons. Since one needs to simulate the shape of each spike, the time step 
in the model must be relatively small, e.g., r = 0.25 ms. Since the model is a 
2-dimensional system of ODEs, without a reset, it cannot exhibit autonomous 
chaotic dynamics or bursting. Adding noise to this, or some other 2-dimensional 
models, allows for stochastic bursting. 



2.2.5 2 component Belousov Zhabotinsky reaction 



Classical Belousov-Zhabotinsky (BZ) reaction is a family of oscillating chemical 
reactions. During these reactions, transition-metal ions catalyze oxidation of 
various, usually organic, reductants by bromic acid in acidic water solution. 
Most BZ reactions are homogeneous. The BZ reaction makes it possible to 
observe development of complex patterns in time and space by naked eye on a 
very convenient human time scale of dozens of seconds and space scale of several 
millimeters. The BZ reaction can generate up to several thousand oscillatory 
cycles in a closed system, which permits studying chemical waves and patterns 
without constant replenishment of reactants [65] . 

Consider the water-in-oil micro-emulsion BZ reaction 35. 34 



d t v = D v Av + — [f z + io (1 - mz)} - — — + — 
£o v + q e 



1 — mz 
1 — mz + e\ 



v — v 



d t z 



D?Az — z + v 



1 — mz 



1 — mz + ei 



(28) 



where v, z are dimensionless concentrations of activator HBr02 and oxidized 
catalyst [Ru(bpy)^\ 3+ respectively; D v and D z are dimensionless diffusion coef- 
ficients of activator and catalyst; /, £o, ej and q are parameters of the standard 
Tyson model [62j : io represents the photoinduced production of inhibitor, and 
m represents the strength of oxidized state of the catalyst with < mz < 1. 
This reaction was shown experimentally and numerically to admit localized spot 
patterns that persist for long time 
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We can rescale the variables in 



as 



z = 1/m — m ' wei, v — m~ 1 v, t = eom 1 t. 

In the new variables, after dropping the hats, we obtain the non-dimensional 
2-component BZ reaction 



dtV = e 2 Av + f(v, w), rdtw — DAw + g(v, w), 



where 



= _[/„ + / lU ,]l-i + 
v + q 



w 



1 + aw 

with the non-dimensional constants given by 

,-V2 



v — v 2 , g(v, w) = 1 



w 



1 + aw 



v, 



a = m 



1/2 / ■ fo 
fi = eim 7 Uo 



1 £ 

e 2 = £ D v m 1/2 , D = D z eimT 1/2 , t = -. 

m £o 

2.3 3— component and multi— component systems 
2.3.1 Oregonator 

The Oregonator model is based on the so-called FKN-mechanism (TT|, which 
provided the first successful explanation of the chemical oscillations that occur 
in the experimental Belousov-Zhabotinsky reaction. It is is composed of five 
coupled elementary chemical stoichiometrics. During the last two decades, the 
Oregonator model has been modified in many ways by inclusion of additional 
chemical reaction steps or by changing the rate constants. If we denote the 
concentration of the species S by [S], then we define: A — [BrO^], H = [H + ], 
X = [HBr0 2 ], Y = [Br"], Z = [Ce 4+ ]. The original Oregonator model was 
described by the following three coupled nonlinear PDEs, 



dtX 
d t Y 
d t Z 



= kiAH 2 Y - k 2 HXY - 2k 3 X 2 + k 4 AHX + D X W 2 X, 
= -k x AH 2 Y - k 2 HXY + k 5 fZ + D Y V 2 r Y, 
= 2k 4 AHX - k b Z + D z VlZ, 



(29) 



where / is a stoichiometric factor [IS], ki(i — 1, 5) are rate constants, while 
Dx, Dy, and Dz are the diffusion constants of the species HBr02, Br~, and 
Ce 4+ respectively (for dilute solutions, the diffusion matrix is diagonal). For 
a thorough discussion of the chemistry on which the Oregonator is based, the 
reader is referred to [52] , 

The Oregonator temporal mass-action dynamics is a well-stirred, homoge- 
neous system of ODEs given by 

X = k\AY - k 2 XY + k a AX - 2k 4 X 2 , Y = -k x AY - k 2 XY + l/2kJBZ, 

Z = 2k 3 AX - k c BZ, 
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which are typically scaled as [52] 

e(dx/dr) = qy~xy+x(l — x), eidy/dr) — —qy—xy+fz, dz/dr = x—z. 

The basic chemistry of the BZ-oscillations involves jumps between high and low 
HBr02 (X) states, which is reflected in the relaxation oscillator nature of the 
Oregonator. This fundamental bistability may be stabilized in a flow reactor 
(CSTR) with reactants and Br~ in the feed stream. Hysteresis between the two 
states is observed both experimentally and in the Oregonator. Quasiperiodicity 
and chaos also are observed in CSTR and can be modeled by the Oregonator 

2.3.2 Multi phase tumor growth equations 

Our last reaction-diffusion system is a general model of multi-phase tumor 
growth, in the form of nonlinear parabolic PDE, as reviewed recently in [57j 

d t t> % = V • (A$i) - V ■ (v,$ 4 ) + Ai(* is Ct) - fiii^u Q) (30) 

(d t = d/dt), where for phase i, $j is the volume fraction (£\ $j = 1), Di is 
the random motility or diffusion, \i(<&i,Ci) is the chemical and phase depen- 
dent production, and [ii(<&i,Ci) is the chemical and phase dependent degrada- 
tion/death, and Vj is the cell velocity defined by the constitutive equation 

Vj = -fA7p, (31) 

where fj, is a positive constant describing the viscous-like properties of tumor 
cells and p is the spheroid internal pressure. 

In particular, the multi-phase equation (|30[) splits into two heat-like mass- 
conservation PDEs [57] . 

S t $ c = S c - V • (v c $ c ), S t $ F = S F - V • (v F $ F ), (32) 

where $ and $ F are the tissue cell/matrix and fluid volume fractions, respec- 
tively, v c and v F are the cell/matrix and the fluid velocities (both defined by 
their constitutive equations of the form of (|3"Tj) ), S c is the rate of production of 
solid phase tumor tissue and S F is the creation/degradation of the fluid phase. 
Conservation of matter in the tissue, <I> + 4> F = 1, implies that V-(v c '$ c ' 
+v F $ F ) = $ c + $ F . The assumption that the tumor may be described by 
two phases only implies that the new cell/matrix phase is formed from the fluid 
phase and vice versa, so that S c + S F = 0. The detailed biochemistry of tumor 
growth can be coupled into the model above through the growth term S , with 
equations added for nutrient diffusion, see [57] and references therein. 

The multi-phase tumor growth model (|30p has been derived from the clas- 
sical transport/mass conservation equations for different chemical species [57) . 

d tUi = P t - V • Ni. (33) 
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Here Ci are the concentrations of the chemical species, subindex a for oxygen, 
b for glucose, c for lactate ion, d for carbon dioxide, e for bicarbonate ion, / 
for chloride ion, and g for hydrogen ion concentration; Pi is the net rate of 
consumption/production of the chemical species both by tumor cells and due 
to the chemical reactions with other species; and Nj is the flux of each of 
the chemical species inside the tumor spheroid, given (in the simplest case of 
uncharged molecules of glucose, O2 and CO2) by Fick's law, 

Nj = -DiVuu 

where Di are (positive) constant diffusion coefficients. In case of charged molecules 
of ionic species, the flux Nj contains also the (negative) gradient of the volume 
fractions <I>i. 

There are three distinct stages to cancer development: avascular, vascular, 
and metastatic - researchers often concentrate their efforts on answering spe- 
cific OUPC-related questions on each of these stages [57]. In particular, as some 
tumor cell lines grown in vitro form spherical aggregates, the relative cheapness 
and ease of in vitro experiments in comparison to animal experiments has made 
3D multicellular tumor spheroids (MTS, see Figure 6 in [57]) very popular in 
vitro model system of avascular tumoral [35]. They are used to study how lo- 
cal micro-environments affect cellular growth/decay, viability, and therapeutic 
response [53] ■ MTS provide, allowing strictly controlled nutritional and me- 
chanical conditions, excellent experimental patterns to test the validity of the 
proposed mathematical models of tumor growth/decay [53] , 

3 Dissipative evolution under the Ricci flow 

In this section we will derive the geometric formalism associated with the 
quadratic Ricci-flow equation (|4]) , as a general framework for all presented bio- 
reaction-diffusion systems. 

3.1 Geometrization Conjecture 

Geometry and topology of smooth surfaces are related by the Gauss-Bonnet 
formula for a closed surface S (see, e.g., [251 \2E\ ) 



13 In vitro cultivation of tumor cells as multicellular tumor spheroids (MTS) has greatly 
contributed to the understanding of the role of the cellular micro-environment in tumor biology 
(for review see |59l 1381 ^ ■ These spherical cell aggregates mimic avascular tumor stages or 
micro-metastases in many aspects and have been studied intensively as an experimental model 
reflecting an in vivo-like micro-milieu with 3D metabolic gradients. With increasing size, most 
MCTS not only exhibit proliferation gradients from the periphery towards the center but they 
also develop a spheroid type-specific nutrient supply pattern, such as radial oxygen partial 
pressure gradients. Similarly, MCTS of a variety of tumor cell lines exhibit a concentric histo- 
morphology, with a necrotic core surrounded by a viable cell rim. The spherical symmetry is an 
important prerequisite for investigating the effect of environmental factors on cell proliferation 
and viability in a 3D environment on a quantitative basis. 




(34) 



17 



where dA is the area element of a metric g on E, K is the Gaussian curva- 
ture, x(E) is the Euler characteristic of E and gen(E) is its genus, or number 
of handles, of E. Every closed surface E admits a metric of constant Gaussian 
curvature K = +1, 0, or —1 and so is uniformized by elliptic, Euclidean, or hy- 
perbolic geometry, which respectively have gen(5 2 ) = (sphere), gen(T 2 ) = 1 
(torus) and gen(E) > 1 (torus with several holes). The integral is a topolog- 
ical invariant of the surface E, always equal to 2 for all topological spheres S 2 
(that is, for all closed surfaces without holes that can be continuously deformed 
from the geometric sphere) and always equal to for the topological torus T 2 
(i.e., for all closed surfaces with one hole or handle). 

Topological framework for the Ricci flow @ is Thurston's Geometrization 
Conjecture |60j . which states that the interior of any compact 3-manifold can 
be split in an essentially unique way by disjoint embedded 2-spheres S 2 and tori 
T 2 into pieces and each piece admits one of 8 geometric structures (including (i) 
the 3-sphere S 3 with constant curvature +1; (ii) the Euclidean 3-space R 3 with 
constant curvature and (iii) the hyperbolic 3-space H 3 with constant curvature 
— 1)0 The geometrization conjecture (which has the Poincare Conjecture as 
a special case) would give us a link between the geometry and topology of 3- 
manifolds, analogous in spirit to the case of 2-surfaces. 

In higher dimensions, the Gaussian curvature K corresponds to the Riemann 
curvature tensor 9\m on a smooth n— manifold M , which is in local coordinates 
on M denoted by its (4, 0)— components Rijki, or its (3, 1)— components R\- k 
(see Appendix, as well as e.g., [25l[28]). The trace (or, contraction) of 9\m, in 
(4, 0)— case using the inverse metric tensor g 1 ^ — (gij)^ 1 , is the Ricci tensor *Hc, 
the contracted curvature tensor, which is in a local coordinate system {x*}™ =1 
defined in an open set U C M, given by 

Rij = tr(fHm) = g kl R m 

(using Einstein's summation convention), while the scalar curvature is now given 
by the second contraction of as 

R = tr(fRc) = g ij Rij. 

In general, the Ricci flow gij(t) is a one-parameter family of Riemannian 
metrics on a compact n— manifold M governed by the equation @, which has 
a unique solution for a short time for an arbitrary smooth metric gij on M 
[IB]. If £Hc > at any local point x = {x 1 } on M, then the Ricci flow ^ 



14 Another five allowed geometric structures are represented by the following examples: (iv) 
the product S 2 X S 1 ; (v) the product H 2 X S 1 of hyperbolic plane and circle; (vi) a left invariant 
Riemannian metric on the special linear group SX(2,R); (vii) a left invariant Riemannian 
metric on the solvable Poincare-Lorentz group £7(1,1), which consists of rigid motions of a 
(1 + 1)— dimensional space-time provided with the flat metric dt 2 — dx 2 ; (viii) a left invariant 
metric on the nilpotent Heisenberg group, consisting of 3 X 3 matrices of the form 
1 * * 

1 * .In each case, the universal covering of the indicated manifold provides a 
1 

canonical model for the corresponding geometry 1431 . 
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contracts the metric gij(t) near x, to the future, while if < 0, then the flow 
([2]) expands g%j(t) near x. The solution metric gij{t) of the Ricci flow equation 
^ shrinks in positive Ricci curvature direction while it expands in the negative 
Ricci curvature direction, because of the minus sign in the front of the Ricci 
tensor Rij. In particular, on a 2-sphere S 2 , any metric of positive Gaussian 
curvature will shrink to a point in finite time. At a general point, there will 
be directions of positive and negative Ricci curvature along which the metric 
will locally contract or expand (see [1]). Also, if a simply-connected compact 
3-manifold M has a Riemannian metric gij with positive Ricci curvature then 
it is diffcomorphic to the 3-sphere S 3 [T5]. More generally speaking, the Ricci 
flow deforms manifolds with positive Ricci curvature to a point which can be 
renormalized to the 3-sphcrc. 

3.2 Reaction diffusion— type evolution of curvatures and 
volumes 

All three Riemannian curvatures (R, 9lc and -Dtra), as well as the associated 
volume forms, evolve during the Ricci flow @. In general, the Ricci-flow evo- 
lution equation ^ for the metric tensor g^ implies the reaction-diffusion-type 
evolution equation for the Riemann curvature tensor 5Hm on an n— manifold M, 

d t mm = AfRm + Q n , (35) 

where Q n is a quadratic expression of the Riemann n— curvatures, corresponding 
to the n— component bio-chemical reaction, while the term A£Km corresponds 
to the n— component diffusion. From the general n— curvature evolution (|35| 
we have two important particular cases o 

1. The evolution equation for the Ricci curvature tensor d\t on a 3-manifold 

M, 

<9 t «Kc = A5Hc + Q 3 , (36) 

where Q3 is a quadratic expression of the Ricci 3-curvatures, corresponding to 
the 3-component bio-chemical reaction, while the term A5Hc corresponds to the 
3-diffusion. 

15 By expanding the maximum principle for tensors, Hamilton proved that Ricci flow g(t) 
given by (J2j preserves the positivity of the Ricci tensor IHc on 3-manifolds (as well as of the 
Riemann curvature tensor IHm in all dimensions); moreover, the eigenvalues of the Ricci tensor 
on 3— manifolds (and of the curvature operator <Km on 4— manifolds) are getting pinched point- 
wisely as the curvature is getting large 1 161 1171 . This observation allowed him to prove the 
convergence results: the evolving metrics (on a compact manifold) of positive Ricci curvature 
in dimension 3 (or positive Riemann curvature in dimension 4) converge, modulo scaling, to 
metrics of constant positive curvature. 

However, without assumptions on curvature, the long time behavior of the metric evolving 
by Ricci flow may be more complicated |51] . In particular, as t approaches some finite time 
T, the curvatures may become arbitrarily large in some region while staying bounded in its 
complement. On the other hand, Hamilton 19 discovered a remarkable property of solutions 
with nonnegative curvature tensor JRm in arbitrary dimension, called the differential Harnack 
inequality, which allows, in particular, to compare the curvatures of the solution of J3J at 
different points and different times. 
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2. The evolution equation for the scalar curvature R, 

d t R = Ai? + 2|<Hc| 2 , (37) 

in which the term 2 1 S>^c | 2 corresponds to the 2-component bio-chemical reaction, 
while the term AR corresponds to the 2-component diffusion. By the maximum 
principle (see subsection 13.4)) . the minimum of the scalar curvature R is non- 
decreasing along the flow g(t), both on M and on its boundary dM (see [ST]). 

Let us now see in detail how various geometric quantities evolve given the 
short-time solution of the Ricci flow equation ^ on an arbitrary n— manifold 
M. For this, we need first to calculate the variation formulas for the Christoffel 
symbols and curvature tensors on M; then the corresponding evolution equa- 
tions will naturally follow (see [Jj>] [7J [S]). If g(s) is a 1-parameter family of 
metrics on M with d s gij = Vij , then the variation of the Christoffel symbols 1^ 
on M is given by 

0.r£ = \g kl (Viity + Vjv u - v,« y ) , (38) 

(where V is the covariant derivative with respect to the Riemannian connection) 
from which follows the evolution of the Christoffel symbols T k j under the Ricci 
flow g(t) on M given by 

d t T k l0 = ~g kl (V;% + VjRa - V,%) . 

From (|38[) we calculate the variation of the Ricci tensor Rij on M as 

d s R tJ = v m (d s r%) - v, (dsTZj ) , (39) 

and the variation of the scalar curvature R on M by 

d s R = -A^ + div(divt;) - (v,£Rc) , (40) 

where V = <? y fjj = tr(u) is the trace of v = (uy). 

If an n— manifold M is oriented, then the volume n— form on M is given, in 
a positively-oriented local coordinate system {x 1 } € 17 C M, by [23] 



^ = det(gij) dx 1 A da; 2 A ... A dx™. (41) 

If d s gij — , then d s dfi = ^Vd/i. The evolution of the volume n— form dfi under 
the Ricci flow g(t) on M is given by the exponential decay/growth relation with 
the scalar curvature R as the (variable) rate parameter, 

d t dfx = -Rdfi, (42) 

which gives an exponential decay for R >— a > (elliptic geometry) and 
exponential growth for R <= a < (hyperbolic geometry) - for any small 
constant a (scalar curvature must be bounded away from zero) . The elementary 
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volume evolution (|4"2"|) implies the integral form of the exponential relation for 
the total n— volume 



vol (.9) = / d M> 
Jm 



as 



d t vo\{g{t)) = - f Rdfi, 



Jm 



which again gives an exponential decay for elliptic R > and exponential growth 
for hyperbolic R < 0. 

This is a crucial point for the tumor suppression by the body: the immune 
system needs to keep the elliptic geometry of the MTS, evolving by (f3"CT|) - by 
all possible means|^f| In the healthy organism this normally happens, because 
the initial MTS started as a spherical shape with R > 0. The immune system 
just needs to keep the MTS in the spherical/elliptic shape and prevent any 
hyperbolic distortions of R < 0. Thus, it will naturally have an exponential 
decay and vanish. 

Since the n— volume is not constant and sometimes we would like to prevent 
the solution from shrinking to an n— point on M (elliptic case) or expanding to 
an infinity (hyperbolic case), we can also consider the normalized Ricci flow on 
M (see Figure [1] as well as ref. [7]): 



is the average scalar curvature on M. We then have the n— volume conservation 
law: 



To study the long-time existence of the normalized Ricci flow (|43|) on an 
arbitrary n— manifold M, it is important to know what kind of curvature con- 
ditions are preserved under the equation. In general, the Ricci flow g(t) on 
M, as defined by the fundamental relation ([2]), tends to preserve some kind of 
positivity of curvatures. For example, positive scalar curvature R (i.e., elliptic 
geometry) is preserved both on M and on its boundary dM in any dimension. 
This follows from applying the maximum principle to the evolution equation 
(|3"T|) for scalar curvature R both on M and on dM . Similarly, positive Ricci 
curvature is preserved under the Ricci flow on a 3-manifold M. This is a special 
feature of dimension 3 and is related to the fact that the Riemann curvature 
tensor may be recovered algebraically from the Ricci tensor and the metric on 
3-manifolds [7]. 

16 As a tumor decay control tool, a monoclonal antibody therapy is usually proposed. Mon- 
oclonal antibodies (mAb) are mono-specific antibodies that are identical because they are 
produced by one type of immune cell that are all clones of a single parent cell. Given (almost) 
any substance, it is possible to create monoclonal antibodies that specifically bind to that sub- 
stance; they can then serve to detect or purify that substance. The invention of monoclonal 
antibodies is generally accredited to Georges Kohler, Cesar Milstein, and Niels Kaj Jerne in 
1975 |36| . who shared the Nobel Prize in Physiology or Medicine in 1984 for the discovery. 
The key idea was to use a line of myeloma cells that had lost their ability to secrete antibodies, 
come up with a technique to fuse these cells with healthy antibody producing B-cells, and be 
able to select for the successfully fused cells. 




where 




(43) 



d t vol(g{tj) = 0. 
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♦ • • 

Figure 1: An example of Ricci flow normalization: unnormalized flow (up) and 
normalized flow (down). 

In particular, we have the following result for 2-surfaces (see [H]): Let 
S = dM be a closed 2-surface, which is a boundary of a compact 3-manifold M. 
Then for any initial 2-metric go on dM, the solution to the normalized Ricci flow 
(|43|) on dM exists for all time. In other words, the normalized Ricci flow in 2D 
always converges. Moreover, (i) if the Euler characteristic of dM is non-positive, 
then the solution metric g(t) on dM converges to a constant curvature metric 
as t — > oo; and (ii) if the scalar curvature R of the initial metric go is positive, 
then the solution metric g(t) on dM converges to a positive constant curvature 
metric as t — > oo. (For surfaces with non-positive Euler characteristic, the proof 
is based primarily on maximum principle estimates for the scalar curvature.) 

Applying to the tumor evolution ([50)1 , the normalized Ricci flow P5|) of the 
MTS will make it completely round with a geometric sphere shell, which is 
ideal for surgical removal. This is our second option for the MTS growth/decay 
control. If we cannot force it to exponential decay, then we must try to normalize 
into a round spherical shell - which is suitable for surgical removal. 

The negative flow of the total n— volume vo\(g(t)) represents the Einstein- 
Hilbert functional (see [44jl7ll4]) 

E(g) = f Rdii = -d t vo\(g(t)). 

If we put d s gij = Vy, we have 

d s E(g) = J \ -AV + div(divv)- (v,fRc) + ^RV^dn = (v^Rgq - dp, 

so, the critical points of E(g) satisfy Einstein's equation \Rgij — R%j in the 
vacuum. The gradient flow of E{g) on an n— manifold M, 

d t9ij =2{VE{g)) ij =Rg ij -2R ij , 

is the Ricci flow ^ plus Rag. Thus, Einstein metrics are the fixed points of 
the normalized Ricci flow! 17 ! 

17 Einstoin metrics on n— manifolds are metrics with constant Ricci curvature. However, 
along the way, the deformation will encounter singularities. The major question, resolved by 
Perelman, was how to find a way to describe all possible singularities. 
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Let A denote the Laplacian acting on functions on a closed n— manifold M, 
which is in local coordinates {x 1 } € U C M given by 

A = //'-'V.V, = g» {dtj - T%d k ) . 

For any smooth function / on M we have [16 , 8 

AVj/ = ViAf+RijVrf and A|V/| 2 = 2|V i V i /| 2 +2i^V i /V i /+2V i /V i A/, 

from which it follows that if we have 

mc > 0, A/ = 0, |V/| = 1, then 
VV/ = and 9tc(V/, V/) = 0. 

Using this Laplacian A, we can write the linear heat equation on M as 
dtu = An, where u is the temperature. In particular, the Laplacian acting on 
functions with respect to git) will be denoted by A g M . If (M, g(t)) is a solution 
to the Ricci flow equation @, then we have 

d t A g(t) = 2R ij V i V j . 

The evolution equation (|37[) for the scalar curvature R under the Ricci flow 
follows from (|40|) . Using equation (|5 1 1) from Appendix, we have: 

div(SHc) = -VR, so that div(div(SKc)) = -AR, 

showing again that the scalar curvature R satisfies a heat-type equation with a 
quadratic nonlinearity both on a 3-manifold M and on its boundary 2-surface 
DM. 

Next we will find the exact form of the evolution equation |[55j) for the Ricci 
tensor 9\c under the Ricci flow git) given by ^ on any 3— manifold M. (Note 
that in higher dimensions, the appropriate formula of huge complexity would 
involve the whole Riemann curvature tensor 9tm.) Given a variation d s gij = Uy, 
from (|39|) we get 

d.Rij = i {Amj + ViVjF - ViCdivt;),- - V^diw),) , 

where A l denotes the so-called Lichnerowicz Laplacian (which depends on Dim, 
see PUIS]). Since 

ViVjR- Vi(div(9tc)) j - Vj(div(SRc))i = 0, 

by (|5T|) (after some algebra) we get that under the Ricci flow @ the evolution 
equation for the Ricci tensor 5Hc on a 3-manifold M is 

dtRij = ARij + 3RRij - QRimRjm + (2|«c| 2 - R 2 ) gij . 

So, just as in case of the evolution ([37]) of the scalar curvature dtR (both on 
a 3-manifold M and on its 2-boundary dM), we get a heat-type evolution 
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equation with a quadratic nonlinearity for dtRij, which means that positive 
Ricci curvature (D\c > 0) of elliptic 3-geometry is preserved under the Ricci 
flow g(t) on M. 

More generally, we have the following result for 3-manifolds (see [H]): Let 
(M, go) be a compact Riemannian 3— manifold with positive Ricci curvature 
9^c. Then there exists a unique solution to the normalized Ricci flow g(t) on 
M with g(0) = go for all time and the metrics g(t) converge exponentially fast 
to a constant positive sectional curvature metric g x on M. In particular, M 
is diffcomorphic to a 3-sphere S 3 . (As a consequence, such a 3— manifold M 
is necessarily diffcomorphic to a quotient of the 3— sphere by a finite group of 
isometries. It follows that given any homotopy 3— sphere, if one can show that 
it admits a metric with positive Ricci curvature, then the Poincare Conjecture 
would follow [7].) In addition, compact and closed 3— manifolds which admit a 
non-singular solution can also be decomposed into geometric pieces [20] . 

From the geometric evolution equations reviewed in this subsection, we see 
that both short-time and long-time geometric solution can always be found for 
2-compo-nent bio-reaction-diffusion equations, as they correspond to evolution 
of the scalar 2-curvature R. Regarding the 3-component bio-reaction-diffusion 
equations, corresponding to evolution of the Ricci 3-curvature 5Hc, we can always 
find the short-time geometric solution, while the long-time solution exists only 
under some additional (compactnes and/or closure) conditions. Finally, in case 
of n— component bio-reaction-diffusion equations, corresponding to evolution of 
the Riemann n— curvature 5Hm, only short-time geometric solution is possible. 

3.3 Dissipative solitons and Ricci breathers 

An important class of bio-reaction-diffusion systems are dissipative solitons 
(DSs), which are stable solitary localized structures that arise in nonlinear 
spatially extended dissipative systems due to mechanisms of self-organization. 
They can be considered as an extension of the classical soliton concept in conser- 
vative systems. Apart from aspects similar to the behavior of classical particles 
like the formation of bound states, DSs exhibit entirely nonclassical behavior - 
e.g., scattering, generation and annihilation - all without the constraints of en- 
ergy or momentum conservation. The excitation of internal degrees of freedom 
may result in a dynamically stabilized intrinsic speed, or periodic oscillations of 
the shape. 

In particular, stationary DSs are generated by production of material in the 
center of the DSs, diffusive transport into the tails and depletion of material 
in the tails. A propagating pulse arises from production in the leading and 
depletion in the trailing end. Among other effects, one finds periodic oscillations 
of DSs, the so-called 'breathing' dissipative solitons [T5] . 

DSs in many different systems show universal particle-like properties. To 
understand and describe the latter, one may try to derive 'particle equations' 
for slowly varying order parameters like position, velocity or amplitude of the 
DSs by adiabatically eliminating all fast variables in the field description. This 
technique is known from linear systems, however mathematical problems arise 
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from the nonlinear models due to a coupling of fast and slow modes [12] • 

Similar to low-dimensional dynamic systems, for supercritical bifurcations 
of stationary DSs one finds characteristic normal forms essentially depending on 
the symmetries of the system; e.g., for a transition from a symmetric stationary 
to an intrinsically propagating DS one finds the Pitchfork normal form for the 
DS-velocity v [5], 

v = (er — 0"q)v — |v| 2 v 

where a represents the bifurcation parameter and ctq the bifurcation point. For 
a bifurcation to a 'breathing' DS, one finds the Hopf normal form |26j 

A= {cr- a )A-\A\ 2 A 

for the amplitude A of the oscillation. Note that the above problems do not 
arise for classical solitons as inverse scattering theory yields complete analytical 
solutions [15j . 

3.3.1 Ricci breathers and Ricci solitons 

Closely related to dissipative solitons are the so-called breathers, solitonic struc- 
tures given by localized periodic solutions of some nonlinear soliton PDEs, in- 
cluding the exactly-solvable sine-Gordon equation^] and the focusing nonlinear 
Schrodinger equation^ 

A metric g%j{t) evolving by the Ricci flow g(t) given by ^ on any 3-manifold 
M is called a Ricci breather, if for some t\ < ti and a > the metrics agij(ti) 
and gij{t2) differ only by a diffeomorphism; the cases a = l,a < l,a > 1 
correspond to steady, shrinking and expanding breathers, respectively. Trivial 
breathers on M, for which the metrics gij{ti) and gijfa) differ only by diffeo- 
morphism and scaling for each pair of t\ and t%, are called Ricci solitons. Thus, 
if one considers Ricci flow as a dynamical system on the space of Riemannian 
metrics modulo diffeomorphism and scaling, then breathers and solitons corre- 
spond to periodic orbits and fixed points respectively. At each time the Ricci 

18 An exact solution u = u(x,t) of the (1+1)D sine— Gordon equation 

d t 2U = d x 2it — sinu, is [l] 

y \ — u> 2 cos(iirf) \ 
uj cosh(v / l — <^ 2 x) J 

which, for uj < 1, is periodic in time t and decays exponentially when moving away from 
x = 0. 

19 The focusing nonlinear Schrodinger equation is the dispersive complex— valued (1+1)D 
PDE p], 

i dtu + d x 2U + \u\ 2 u = 0, 
with a breather solution of the form: 

(2b 2 cosh(0) + 2 i b yj2 - b' 2 sinh(6») \ , 2 . 2 / 

«= \- r-^ = --1 Lexpia 2 t with 6 = a 2 b V 2 - b 2 t, 

\ 2 cosh(6») - V2 V2 - b' 2 cos(a b x) J 

which gives breathers periodic in space x and approaching the uniform value a when moving 
away from the focus time t = 0. 



u = 4 arctan 
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soliton metric satisfies on M an equation of the form [5T] 

Rij + C9ij +Vib + Vjbi = 0, 

where c is a number and hi is a 1-form; in particular, when = ^Vja for some 
function a on M, we get a gradient Ricci soliton. An important example of a 
gradient shrinking soliton is the Gaussian soliton, for which the metric gij is 
just the Euclidean metric on R 3 , c = 1 and a — — \x\ 2 /2. 

3.4 Smoothing/Averaging heat equation and Ricci entropy 

Given a C 2 function u : M — > R on a Riemannian 3— manifold M, its Laplacian 
is defined in local coordinates {x 1 } S £/ C M to be 

Am = tr g (V 2 u) = //'-' V, V . //. 

where is the covariant derivative (Levi-Civita connection, see Appendix). 
We say that a C 2 function u : M x [0, T) — ► R, where T S (0, oo], is a solution 
to the heat equation if ([3]) holds. One of the most important properties satisfied 
by the heat equation is the maximum principle, which says that for any smooth 
solution to the heat equation, whatever point-wise bounds hold at t = also 
hold for t > 7;. More precisely, we can state: Let u : M X [0, T) — * R be 
a C 2 solution to the heat equation |3|) on a complete Riemannian 3— manifold 
M. If Ci < u(x,0) < C 2 for all x e M, for some constants Ci,C 2 € R, then 
Ci < u(x,t) < C 2 for all x e M and £ S [0,T). This property exhibits the 
averaging behavior of the heat equation Q on M. 

Now, consider Perelman's entropy functional |51j on a 3-manifold Ail 2 "! 

^ = / (i2+|V/| 2 )e-^^ (44) 

J A/ 

for a Riemannian metric gij and a (temperature-like) scalar function / (which 
satisfies the backward heat equation) on a closed 3-manifold M, where d/i is 
the volume 3-form (j4Tj) . During the Ricci flow ([2]), T evolves on M as 

d t T = l\ \Rij + ViVjf\ 2 e-fdfi. (45) 



"Note that in the related context of Riemannian gravitation theory, the so— called gravi- 
tational entropy is embedded in the Weyl curvature (4,0)— tensor 20, which is the traceless 
component of the Riemann curvature tensor 5Hm (i.e., iKm with the Ricci tensor IHc removed), 

2H= Vtoi-fiRijgij), 

where f(Rijgij) is a certain linear function of Rij and gy. According to Penrose's Weyl 
curvature hypothesis, the entire history of a closed universe starts from a uniform low-entropy 
Big Bang with zero Weyl curvature tensor of the cosmological gravitational field and ends with 
a high— entropy Big Crunch, representing the congealing of may black holes, with Weyl tensor 
approaching infinity (see 51) ST]). 
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Now, define X(gij) — inf F{gij, /), where infimum is taken over all smooth /, 
satisfying 

/ e~ f diJ, = l. (46) 
Jm 

X(gij) is the lowest eigenvalue of the operator — 4 A + R. Then the entropy 
evolution formula (|45[) implies that \{gij (t)) is nondecreasing in t, and moreover, 
if X(ti) = A(<2)) then for t £ [tijia] we have Rij + ViVjf = for / which 
minimizes T on M [51j . Thus a steady breather on M is necessarily a steady 
soliton. 

If we define the conjugate heat operator on M as 

□* = -d/dt-A + R 

then we have the conjugate heat eauatioi^\ [51] 

□*u = 0. (48) 

The entropy functional (|44[) is nondecreasing under the coupled Ricci-diffusion 
flow on M (see [BSEO]) 

dtgij = -2i2y, 9 t u = -Ati + -u - ■ L ^- L , (49) 

where the second equation ensures J M u 2 d/i = 1, to be preserved by the Ricci 

flow g{t) on M. If we define u = e~i , then the right-hand equation in (|49p is 
equivalent to the generic scalar-field /—evolution equation on M, 

d t f=-Af-R+\Vf\ 2 , 

21 In [51) Perelman stated a differential Li-Yau-Hamilton (LYH) type inequality |22J for the 
fundamental solution u = u(x, t) of the conjugate heat equation (148 I t on a closed n— manifold 
M evolving by the Ricci flow (f2j) . Let p £ M and 

u = (4ttt)~ ~2e~f 

be the fundamental solution of the conjugate heat equation in M X (0, T), 

□*u = 0, or dtu + An = Ru, 

where r = T — t and R = R(-, i) is the scalar curvature of M with respect to the metric g(t) 
with liirii y<T u = S p (in the distribution sense), where <5 P is the delta-mass at p. Let 

v = [t(2A/ - I V/| 2 + R) + f- n]u, 

where r = T — t. Then we have a differential LYH-type inequality 

v(x,t)<0 in Mx(0,T). (47) 

This result was used by Perelman to give a proof of the pseudolocality theorem |51l which 
roughly said that almost Euclidean regions of large curvature in closed manifold with metric 
evolving by Ricci flow g(i) given by J2} remain localized. 

In particular, let (M, g(t)), < t < T, dM 7^ <j>, be a compact 3— manifold with metric g(t) 
evolving by the Ricci flow g(t) given by J2J such that the second fundamental form of the 
surface dM with respect to the unit outward normal d/dv of dM is uniformly bounded below 
on dM X [0, T] . A global Li-Yau gradient estimate 1391 for the solution of the generalized 
conjugate heat equation was proved in 1221 (using a a variation of the method of P. Li and 
S.T. Yau, |39p on such a manifold with Neumann boundary condition. 
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which instead preserves PS]) . 

The coupled Ricci-diffusion flow ([49]) , or cquivalently, the dual system 

d m A-.;'/,, • Q.Ag. <),/). d t f = -Af-R+\Vf\ 2 , (50) 

is our global decay model for a general n— dimensional bio-reaction-diffusion 
process, including both geometric and bio-chcmical multi-phase evolution. 

The sole hypothesis of this paper is that any kind of reaction-diffusion pro- 
cesses in biology, chemistry and physics is subsumed by the geometric-diffusion 
system (|49|). or the dual system ([50]) . 

4 Appendix: Riemann and Ricci curvatures on 
a smooth n— manifold 

Recall that proper differentiation of vector and tensor fields on a smooth Rie- 
mannian n— manifold is performed using the Levi-Civita covariant derivative 
(see, e.g., [25] HE])- Formally, let M be a Riemannian n— manifold with the 
tangent bundle TM and a local coordinate system {x l }™ =1 defined in an open 
set U C M. The covariant derivative operator, V x ■ C°°{TM) -> C°°(TM), is 
the unique linear map such that for any vector fields X, Y, Z, constant c, and 
function / the following properties are valid: 

Vx +C y = V x +cVy, 
S7 x {Y + fZ) = V x Y+(Xf)Z + fV x Z, with 
Vjf Y — VyX = [X,Y], (torsion free property) 

where [X, Y] is the Lie bracket of X and Y (see, e.g., [23]). In local coordinates, 
the metric g is defined for any orthonormal basis (di = d x i ) in U C M by 

9ij = f}{9i,dj) = 5ij, d k gij = 0. 

Then the affine Levi-Civita connection is defined on M by 

V„ dj = T% <h, , where T% = X -g kl (d i9jl + d j9a - d m ) 

are the (second-order) Christoffel symbols. 

Now, using the covariant derivative operator Vx we can define the Riemann 
curvature (3, 1)— tensor fHm by (see, e.g., [25l |28] 1 

mm(X, Y)Z = VxVyZ - VyVxZ - V [X . Y ]Z. 

SHm measures the curvature of the manifold by expressing how noncommutative 
covariant differentiation is. The (3, 1)— components R\j k of D\m are defined in 
UcMby 

y\m(di,dj)dk = R\j k di, which expands fsee [44] ) as 

Til pt yl fi yl i pm-pi pmri! 

n ijk — °i l jk ~ °J l ik 1 jk L im ~ 1 ik l jm- 
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Also, the Riemann (4,0)— tensor Rijki — gimR^k 1& defined as the g— based 
inner product on M, 

Rijkl = (di,dj) d k ,di) . 

The first and second Bianchi identities for the Riemann (4,0)— tensor R^ki 
hold, 

Rijkl + Rjkil + Rkijl = 0, ViRjklm + V ' jRkilm + V ' kRijlrn = 0, 

while the twice contracted second Bianchi identity reads 

2VjR VJ = ViJZ. (51) 
The (0, 2) Ricci tensor 9lc is the trace of the Riemann (3, 1)— tensor Dtm, 
<ftc(F, Z)+tr(X -> $Km(X, F)Z), so that SHc(X, F) = g(9lm(9,,X)0j, F), 
Its components ijjfc = 5Hc dk) are given in [/ C M by the contraction [44] 
i?jfc = -R^fe, or, in terms of Christoffel symbols, 

p n *p? o -pi _j_ pi pm pi nm 

-fi-jfc — u i L jk u k L ji ' L mi L jk L mk L ji- 

Being a symmetric second-order tensor, has n + 12 independent components 
on an n— manifold M. In particular, on a 3-manifold, it has 6 components, and 
on a 2-surface it has only the following 3 components: 

Rn = g 22 R2ii2, R12 = g 12 R2i2ii R22 — g R1221, 

which arc all proportional to the corresponding coordinates of the metric tensor, 

Rll _ Rl2 _ R22 _ R\2\2 

ffii 9i2 922 det(g)' 

Finally, the scalar curvature R is the trace of the Ricci tensor 5Hc, given in 
UcMby. R = g i iR l] . 
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